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Abstract 

Because of physical processes ranging from microscopic particle collisions to macroscopic 
hydrodynamic fluctuations, any plasma in thermal equilibrium emits gravitational waves. 
For the largest wavelengths the emission rate is proportional to the shear viscosity of the 
plasma. In the Standard Model at T > 160 GeV, the shear viscosity is dominated by the 
most weakly interacting particles, right-handed leptons, and is relatively large. We estimate 
the order of magnitude of the corresponding spectrum of gravitational waves. Even though 
at small frequencies (corresponding to the sub-Hz range relevant for planned observatories 
such as eLISA) this background is tiny compared with that from non-equilibrium sources, 
the total energy carried by the high-frequency part of the spectrum is non-negligible if the 
production continues for a long time. We suggest that this may constrain (weakly) the highest 
temperature of the radiation epoch. Observing the high-frequency part directly sets a very 
ambitious goal for future generations of GHz-range detectors. 
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1. Introduction 


Gravitational waves offer a possible way to observe phenomena taking place very early in 
the history of the universe. Famously, long-wavelength waves are produced during a period 
of inflation (cf. e.g. refs. [1, 2]). However at higher frequencies gravitational waves can also 
probe post-inflationary non-equilibrium phenomena, such as preheating [3]—[9], topological 
defects [10]—[18], bubble dynamics related to a first-order phase transition [19]—[24], or noisy 
turbulent motion [25]—[29]. Recent numerical simulations start to account for both bubble 
dynamics and the subsequent motion [30, 31]. For a review concerning post-inflationary 
sources and the associated observational prospects, see ref. [32], 

It is well-known that gravitational waves are being produced in thermal equilibrium as 
well (for an example, see ref. [33]). In thermal equilibrium particles scatter on each other, 
which implies the presences of forces and accelerations. However, for a physical momentum 
k > T, the thermal production rate is suppressed by e fc / T , because the energy carried 
away by the graviton must be extracted from thermal fluctuations. Since typical particle 
momenta are ~ 3 T and scatterings are proportional to the coupling strengths responsible for 
the interactions, it may be assumed that the rate is suppressed by ~ aT 3 e~ k ^ T /m P1 , where 
a is a fine-structure constant and m pl is the Planck mass. In a weakly-coupled system such 
as the Standard Model, this rate is small. 

On the other hand, a thermal system experiences also long-wavelength fluctuations not 
associated with single particle states. At the smallest momenta l;<T these can be called 
hydrodynamic fluctuations [34], We are not aware of an estimate of a corresponding contri¬ 
bution from equilibrium Standard Model physics to the gravitational wave background, and 
one purpose of the present note is to provide for one. (Very similar physics, although boosted 
by a conjectured turbulent cascade, has recently been discussed in ref. [29].) In contrast to 
the non-equilibrium phenomena mentioned above, the “advantage” of a thermal contribution 
is that it is guaranteed to be present. Another purpose of our study is to roughly estimate 
the production rate at k > T, and to motivate the need for a complete computation. 

In order to be more specific, consider the closely analogous case of the production rate of 
photons from a plasma which is neutral but has electrically charged constituents. Even though 
the expectation values of the electromagnetic charge density and current vanish, (n em ) = 0, 
(J em ) = 0, thermal fluctuations do induce charge fluctuations which have a non-zero root 
mean square value: 
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where V denotes the volume, x e m a susceptibility, and (...) a thermal expectation value. The 
susceptibility is non-zero even without interactions, for instance for a plasma of free massless 
Dirac fermions representing electrons and positrons it reads x e m = e 2 T 2 /3, e 2 = 47ra em . 
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Because of diffusion, the charge fluctuations induce electromagnetic currents, and currents 
in turn source photons. Currents can also directly originate from fluctuations. Assuming 
that the photons produced do not equilibrate as fast as the plasma, which is the case for 
instance for the plasma generated in heavy ion collision experiments, the thermal average of 
their production rate can be evaluated. A text-book computation shows that the rate per 
unit volume can be expressed as [35, 36] 
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where /C = (k, k), k = |k|; X = (t, x); 1C - X = kt — k • x; and denote polarization vectors. 

12 

For k = ke%, the polarization sum only couples to the transverse components J ei h. 

For small icT, operator ordering plays no role in eq. (1.2), and the fluctuations are also 
uncorrelated in space and time [34]. Their amplitude is related to diffusion or, equivalently, 
to conductivity (er). This yields finally 


dl\(k) k<alT jTa _ a em T 2 

d 3 k (2ir) 3 k (2n) 3 k a 2 ln(l/a s ) 

where we inserted the parametric form of the conductivity of a QCD plasma [37, 38]. For large 
k>3T, in contrast, the rate originates from particle scatterings rather than hydrodynamic 
fluctuations, and has the parametric form [39, 40] 

dF 7 (k) k>3T a em a s ln(l/ a s )T 2 e~ k ^ T 

“SOT ~ -(2 Wh -’ (L4) 

In the following we show that results analogous to eqs. (1.3) and (1.4) apply to gravitational 
waves, just with the replacements a em —> T 2 /m|j and a s —> a. 

Our presentation is organized as follows. After deriving an expression for the gravitational 
wave production rate in sec. 2, we analyze the structure of the energy-momentum tensor 
correlator for k <C T in sec. 3. The quantity parametrizing this structure, the shear viscosity, 
is briefly discussed in sec. 4. In sec. 5 we turn to the other case k > 3 T and compute the 
logarithmically enhanced terms in this regime. The results are embedded in a cosmological 
background in sec. 6 and compared with a well-studied non-equilibrium source in sec. 7. 
Section 8 offers some conclusions and an outlook. 


2. Production rate of gravitational waves from thermal equilibrium 

As a first ingredient, we consider the rate at which energy density is emitted in gravitational 
waves. The derivation can be carried out in two different ways: by treating gravitons as 
quantized particles, or through a purely classical analysis. We start with the first method, 
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leading to a result analogous to eq. (1.2). We work first in Minkowskian spacetime, adding 
cosmological expansion in sec. 6. 

The linearized equation of motion for the metric perturbation h- in the traceless transverse 
gauge reads 

^-V 2 ^ = 16ttG^ t , (2.1) 

where G = 1/m^. The right-hand side of this equation plays the role of the electromagnetic 
current in the photon case. The classical energy associated with gravitational waves reads 

Egw = 32ttG J xev ’ (2 ' 2) 

where V is a volume. It is well-known that the corresponding energy density cannot be 
localized. However, if we express a free hj? as a usual linear combination of forward and 
backward-propagating plane waves, and omit fast oscillations exp(±2 iut), then eq. (2.2) can 
be re-interpreted as a Hamiltonian with a familiar canonical form: 

H = ((E gw )) = {[g T (*,x)] 2 + |V^ T (*,x)| 2 } . (2.3) 


Here ((...)) denotes an average over an oscillation period. From eq. (2.3) canonically normal¬ 
ized fields can be identified as h 1 T = hj ?/V 32ttG. According to eq. (2.1) they are sourced 
as d^hJJ — V 2 /iT r = V 8nGTJ-'. We can now directly overtake eq. (1.2) for the production 
rate of gravitons, by replacing the polarization vectors accordingly. Subsequently, weighting 
the production rate by the energy carried by individual quanta, we obtain 
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The sum over the polarization vectors yields 
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where P V] = 5- — kikj/ k 2 . Choosing k = ke^ and rotating subsequently the diagonal 
correlator (1 (T 11 — T 22 ) (T 11 — T 22 )) into the non-diagonal one, we obtain 

a nS* = S/* e ' MM,<Tn(0)r,aW> - (2 ' 6) 

In order to be convinced that eq. (2.4) is correct, let us repeat the analysis on a purely 
classical level. Fourier transforming eq. (2.1) with respect to spatial coordinates and denoting 
the retarded Green’s function related to the time evolution by A (t, k), its time derivative reads 
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A (t — t', k ) = Q(t — t') cos (k(t — t')) . Dividing by volume, the energy density corresponding 
to eq. (2.2) can then be expressed as 


Pgw = ^- I I <W dt' , cos(k(t-t'))cos(k(t-t , '))(T^ T (t\'k)T^ T (t' , ,-k)^, (2.7) 


k J — oo J — oo 


where f k = f d 3 k/(27r) 3 . Following now a standard argument, let us assume that the sources 
switch off before the observation time t. Then the upper bounds of the time integrals can be 
treated as independent of t, and we can average over fast oscillations within the integrand: 
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Taking a time derivative and assuming that 4> is a function of the time difference 1 yields 


m = 
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Going finally back to configuration space and making use of translational invariance in spatial 
and temporal directions we get 


dp, 


GW 


(«-"■*> (ip’-Vt.x), 7^(0,0)}) . (2.10) 


dt d 3 k (2 

Given that eq. (2.5) defines a projection operator to the TT modes and that in the classical 
limit operator ordering plays no role, eq. (2.10) indeed agrees with eq. (2.4) for k <C T. 


3. Correlation function in the tensor channel 

Having obtained eq. (2.6), the next task is to determine the shape of the energy-momentum 
tensor correlator in momentum space. Here we do this for small light-like four-momenta 
(w, k < a 2 T), returning to the regime k>3T in sec. 5. 

1 This can be justified, for instance, if the typical time differences are short (say, reflecting physics much 
within the horizon) compared with the observation time scale (say, the Hubble time). 
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Consider hydrodynamic fluctuations associated with a local flow velocity v l around an 
equilibrium state at a temperature T. To first order in gradients and in v \ 2 the energy- 
momentum tensor has the form 


T 01 = (e + p)v\ (3.1) 

T ij = (p - CV • v) S ij - ri(d iV j + djj - \ <WV • v) , (3.2) 

where e,p,(, V are the energy density, pressure, bulk viscosity, and shear viscosity, respec¬ 
tively. The equation for energy-momentum conservation asserts that d 0 T 0j + d{T l i = 0 Vj € 
{1, 2,3}. Let us consider a plane wave perturbation with a momentum vector k = ke 3 . Then 
the equations of motion for the transverse velocity components (v_l • k = 0) decouple from 
the equations relating u 3 and d 3 p. The resulting system is immediately integrated to obtain 

v ± (f, k) = v_|_(0, k) e -^' 2 t /( e +P) . (3.3) 


We now consider the 2-point correlator 

(±{r w (i,k),T°*(0,-k)}} , (3.4) 


where the operator ordering is only relevant in the quantum theory. This correlator is sym¬ 
metric in t —>• —t and has a classical limit. Therefore equations (3.1) and (3.3) lead to a 
hydrodynamic prediction for the transverse components (i',f € { 1 , 2 }), 
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where we returned to configuration space for the equal-time correlator. 

Let us take k to be very small, and look for the leading term in this limit. In the equal¬ 
time correlator we can send k —* 0. Then it equals the susceptibility related to the total 
momentum in the ^-direction: 

[ (T 0i '(0,x)T 0 ^(0,0)) =1/ (T 0i '(0,x)T 0 ^(0,y)) = Xp . (3.6) 

4xev x /TV Jx,yev N ' 


Even though the average momentum is zero, its susceptibility is non-zero, in analogy with 
eq. ( 1 . 1 ): 

X P = T(e+p). (3.7) 

Despite including an integral over operator correlations at short distances, this exact equation 
is ultraviolet finite just like eq. (1.1) (for a rigorous discussion see ref. [41]). 

2 Second order terms such as (e + p)vV are omitted. 
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Now, a Ward identity related to energy-momentum conservation asserts that 


a; 2 (^{T°*'(a;,k),r 0 ^(-a;,-k)}) = k 2 {T 3 i \uj, k), T**'{-u, -k)}) . (3.8) 

Therefore, eqs. (3.5)-(3.7) can be re-expressed as 




Lj,k ^ a 2 T 


2 riTup-fr'i' 
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)? 2 fc 4 

(e+pF 


(3.9) 


Taking lirn^^Q limfc^o from here yields the well-known Kubo formula for ij. 

Of interest to us is not the correlator of eq. (3.9) (known as the “shear channel”) but the 
corresponding correlator for the spatial components transverse to k (known as the “tensor 
channel”). It can be argued, however, that its functional form is closely related to that in 
eq. (3.9). The tensor components also experience hydrodynamical fluctuations; but, in our 
coordinate system with k = ke 3 , these are related to the velocity gradients d\v 2 or 82 V 1 , 
cf. eq. (3.2). These components are decoupled from the equations of motion following from 
energy-momentum conservation. Therefore, they are not represented by smooth differen¬ 
tiable functions responsible for the transfer of hydrodynamic information from one point and 
time to another; rather, nearby points are uncorrelated, as is the case for generic thermal 
fluctuations [34]: 


{^{ T Vf( t i< x i)>2fc'7'(*2,x 2 )}) = -* 2 )£ ( 3 ) ( xi - x 2 ) . (3.10) 


Consequently a Fourier transform like in eq. (3.9) is independent of co, k. Putting finally 
oj = k and sending k —» 0 so that the distinction between spatial directions disappears, we 
can fix the coefficient through a comparison with eq. (3.9): 

j™ j x eik ^ z \^{T 12 (X),T 12 (0)}) =2 V T . (3.11) 


This is the main result that is needed below. We remark that, apart from the physical 
arguments discussed above, the same expression can be derived more formally by a linear 
response analysis related to a metric perturbation (cf. e.g. ref. [42]). 


4. Estimate of shear viscosity at T > 160 GeV 

Shear viscosity (rf) is a macroscopic property of a plasma that originates from the microscopic 
collisions that its constituents are undergoing. It is inversely proportional to a scattering 
cross section and therefore large for a plasma in which there are some weakly interacting 
particles. In the Standard Model above the electroweak crossover, right-handed leptons 
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k ~ 3 T 


X V ' / ' / Vvw $$$& 


p ~ 3T 


Figure 1: Processes leading to a logarithmically enhanced graviton production rate. Wiggly lines 
denote gauge bosons; arrowed lines fermions; dashed lines scalars; and a double line a graviton. By 
k,p ~ 3T we denote typical momenta of the scattering particles, whereas the filled blob indicates that 
the vertical rung carries a soft spacelike momentum transfer (t ~ — ~ — g 2 T 2 , where ■ k = 0) 

so that the gauge boson needs to be Hard Thermal Loop resummed. 


are the most weakly interacting degrees of freedom, changing their momenta only through 
reactions mediated by hypercharge gauge fields. 

Omitting for the moment all particle species which equilibrate faster than right-handed 
leptons, the shear viscosity can be extracted from refs. [37, 38]: 


16T 3 

11 ~ gjln(5T/m D1 ) ’ 


(4.1) 


where m D1 = y/11/6 g\T is the Debye mass related to the hypercharge gauge field. Inserting 
gi ~ 0.36 for the gauge coupling we obtain 


T] ~ 400 T 3 . 


(4.2) 


We use this value for order-of-magnitude estimates below. 

If we increase the temperature above 160 GeV, the hypercharge coupling g\ grows and 
the weak and strong couplings # 2,53 decrease. Presumably, the top Yukawa coupling h t and 
the Higgs self-coupling A are also of a similar magnitude. In this situation the analysis of 
refs. [37, 38] should be generalized to include a scalar field and a more complicated set of 
reactions. Even though conceptually straightforward, implementing and solving numerically 
the corresponding set of rate equations is a formidable task and beyond the scope of the 
present investigation. We note, however, that the shear viscosity is likely to decrease with 
increasing g \, so that eq. (4.2) should represent the most “optimistic” estimate from the point 
of view of detecting a thermally emitted low-frequency gravitational wave background. 


5. Leading-logarithmic production rate at large momentum 

Before turning to numerical estimates we wish to complete the qualitative picture concerning 
the thermal graviton production rate by considering the case of “hard momenta”, k ~ 3 T. A 
full computation of the rate in this regime represents a complicated task, similar to the full 
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computation of the shear viscosity when all couplings are of the same order of magnitude, 
and is postponed to future work. In contrast to the shear viscosity, for hard momenta the 
result is dominated by the largest couplings, in particular the strong gauge coupling. If we 
restrict to logarithmically enhanced terms (cf. eq. (1.4)) then it can be shown that only the 
gauge couplings (< 71,52 > 53 ) contribute at leading order. 

An elegant way to determine the logarithmically enhanced terms has been discussed in 
ref. [43], sec. 4.2. Scatterings experienced by soft space-like gauge bosons, the vertical rung 
in fig. 1, correspond to Landau damping, and can be represented within the Hard Thermal 
Loop (HTL) effective theory [44, 45]. Computing the 2-point correlator of t 12 within the 
HTL theory and noting that only one of the gauge bosons attaching to the graviton vertex 
can be soft at a time, yields (for q <C k ~ 3T) 
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(5.1) 


where f B is the Bose distribution; gy = q • k /k; m 2 is a Debye mass squared; A indicates 
that this treatment only applies to soft modes q± <C 3T; and p T/E are spectral functions 
corresponding to the “transverse” and “electric” polarizations, respectively. In the last step 
we made use of a sum rule for the HTL-resummed gluon propagator that has been derived 
in refs. [46, 47]. 

The integral in eq. (5.1) happens to be identical to that appearing in the context of the jet 
quenching parameter q in QCD [47]. Carrying it out and inserting the result into eq. (2.6), 
we obtain 


d/^GW 
d t d In k 
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where d\ = 1, = 3, d 3 = 8 ; m Di is the Debye mass corresponding to the gauge group 
U(l), SU(2) or SU(3), respectively; g 2 € {g 2 , g 2 , g%, h 2 }', and the ultraviolet scale within 
the logarithm has been (arbitrarily) taken over from eq. (4.1). The Debye masses read 
rn 2 n = llgfT 2 /6 , m^ 2 = llg 2 T 2 /6, and m 2 M = 2^|T 2 . Because of the largest multiplicity, 
the result is dominated by the QCD contribution. We note that a similar computation for 
f-channel fermion or Higgs exchange does not lead to logarithmic enhancement. 


6. Embedding the result in cosmology 


Combining eqs. (2.6), (3.11) and (5.2), we get 
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Figure 2: The result from eq. (6.8), multiplied by TO pl /T max , as a function of the present-day 
frequency. The maximum of the power lies in the range k ~ T max at T = T max , and at k ~ To at 
T = Tq. The hydrodynamic and leading-log results correspond to the two limits shown in eq. (6.2), 
with the band originating from varying ?) = 100...400 in the hydrodynamic prediction and from varying 
the constant 0(1) within the range 0...10 in the leading-logarithmic result. The couplings were fixed 
at a scale g = ttT with T ~ 10 6 GeV: g\ ~ 0.13, g\ ss 0.40, g\ ss 1.0. For obtaining the current day 
energy fraction the result needs to be multiplied by fi , ~5x 10” 5 . The eLISA sensitivity peaks at 
/ ~ 10 -2 ...10 -3 Hz. 


This applies in a local Minkowskian frame. The function 0, 
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is quantitatively correct at k < a 2 T whereas at k > 3 T it only represents the qualitative struc¬ 
ture (in particular the coefficient “5” inside the logarithm is but a convention, and there could 
be substantial non-logarithmic contributions from h 2 or from 0(g)-suppressed effects like in 
the case of the jet quenching parameter q [47]). We would now like to re-express the result 
in an expanding cosmological background, and subsequently obtain numerical estimates. As 
our reference temperature we take that corresponding to the electroweak crossover in the 
Standard Model, Tq = 160 GeV (cf. e.g. refs. [48, 49]). 
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The basic equations needed from cosmology are (for a fiat spatial geometry) 


18ne 1 
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V 3 m pi 
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(6.3) 


where H is the Hubble rate, a is the scale factor, and s(T ) is the entropy density. Combining 
the two equations in eq. (6.3), the relation between time and temperature can be expressed 
as 

^ = -TH(T)3c 2 (T) , (6.4) 

where c s is the speed of sound, c 2 (T) = p'(T)/e'(T). The energy density carried by grav¬ 
itational waves is of the form p GW (t ) = f k k where / is a phase space distribution. 

Making use of the known evolution equation for / in an expanding background, the energy 
density can be seen to evolve as 


(d t + 4H)p GW (t) = f R(T,k) , (6.5) 

J k 

where R(T,k) = 32TnjT(j)(k/T) /m^ in the notation of eq. (6.1). Given that (dt + 3 H)s = 0, 
the factor 4 H can be taken care of by normalizing p GW by s 4 / 3 . Subsequently the equation 
can be integrated, by assuming that at an initial time f min (corresponding to a maximal 
temperature T max ) there were no (thermally produced) gravitational waves present: 


Pgw(*o) = ft* f R(T, k) = /- rmax dT f R(T,k ) 

s 4 / 3 (io) Jt . k sW(t) J To J^TH(T)3c 2 s (T)s^(T)- 


( 6 . 6 ) 


Taking into account that momenta redshift as k(t ) = koa(to)/a(t) and expressing the mo¬ 
mentum space integrals in terms of ko finally yields 
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(6.7) 


where we also inserted H from eq. (6.3). If we approximate c 2 ~ 1/3; assume all thermody¬ 
namic functions to scale with their dimension (s = sT 3 , r] = ?)T 3 , e = eT 4 , with s, fj , e roughly 
constants at T > To); and consider T max S> Tq, then 


^gw(^’o) - 


24fj T max fc| 

V67r 3 e 3 to pi 



( 6 . 8 ) 


This result is plotted in fig. 2, after a redshift of ko/To to a current-day frequency. 

For a given mode To, production starts at a maximal temperature when the argument of 4> 
is of order unity; since the entropy density roughly scales with T 3 for T > Tq, this poses no 
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particular constraint if we restrict to ko < To. The horizon radius of a given period redshifts, 
and we are interested in causal physics taking place within the horizon. For instance, a 
temperature leading to a horizon radius comparable to the planned eLISA [50] arm length, 
~ 10 6 km, is T max ~ 10 6 GeV. Let us take this as an example. Inserting i) ~ 400, e ~ 35, 
T max ~ 10 6 GeV into eq. (6.8), we thus get 

!*(!). 

The Hubble radius ( H _1 ) of the electroweak epoch (T = Tq) corresponds to ~ 10 10 km today, 
so if we consider wavelengths extending up to ~ 10 6 km, then ko ~ aH with a ~ 2ir x 10 4 . 
In this situation Lq/Tq can be estimated as 


k 0 a.H{To) 

- r^j - 

To To 


Sire To 
3 m 


■ a 


2 x 10 -16 a . 


pi 


( 6 . 10 ) 


Inserting into eq. (6.9) we find a very small energy fraction. However the fraction is larger if 
we consider the total amount of energy in gravitational waves, which originates dominantly 
from k 0 ~ ?b; this energy is constrained to be below that corresponding to one equilibrated 
relativistic degree of freedom [51, 52], f ]n kg Q GW ^ 1/100. We return to this consideration 
around eq. (8.2), but first compare the infrared part with non-equilibrium processes. 


7. Order-of-magnitude comparison with a non-equilibrium source 

In order to get a qualified impression about the magnitude of the thermal background, con¬ 
sider the well-studied case of a first-order phase transition at the electroweak epoch. As 
before 9 GW denotes the ratio of the energy densities of gravitational waves and radiation at 
To ~ 160 GeV. Today, the energy density in radiation corresponds to 9 rad ~ 5 x 10~ 5 , and 
the projected eLISA sensitivity is 9 eLISA ~ 10 -11 . So, in order to be detectable, we hope to 
find a signal in the range H GW ~ 2 x 10~ 7 at the electroweak epoch. 

Apart from the overall magnitude, an important feature of any observatory is that its 
sensitivity peaks in a certain frequency range. For eLISA, this is / (10 3 ... 10 2 ) Hz. This 

corresponds to a distance scale l B ~ (10~ 3 ... 10 -2 ) (t K in terms of the horizon radius of the 
electroweak epoch, where we have defined = if -1 . 

Now, the overall signal from a first-order transition is traditionally argued to be of the 
form [19]—[24] 

(^) ’ t 7 ' 1 ) 

where v w is a bubble wall velocity; n > 0; L is the latent heat released in the transition; 
k < 1 parametrizes the efficiency at which energy is converted into gravitational waves [22]; 
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and £ B is the typical bubble separation. The spectrum peaks at momenta corresponding to 
the bubble separation scale, k max ~ 2ir/£ B . For k < k max , a behaviour ~ k 3 has been found, 
whereas at k > A: max the spectrum falls off, perhaps as 1/k or 1/k 2 [24], A recent numerical 
study shows that the complicated dynamics following bubble collisions continues for a long 
time and thereby boosts eq. (7.1) by a factor ~ £ W _/£ B , with the price of a more rapid decay 
of the power spectrum at large k [30]. 

If we optimistically insert £ B ~ 0.0U? H into eq. (7.1), it still remains a challenge to get 
fl GW ~ 10 - ' out. Basically, a very large latent heat L/e > 0.01 would be needed. This 
requires a drastic modification of the Higgs sector, which normally carries just a handful of 
degrees of freedom in comparison with ~ 100 contributing to e. However, with somewhat less 
drastic assumptions, and including a boost by ~ £ H /T B from ref. [30], numbers like H GW ~ 
10~ 10 could be obtained, which is also of interest for future generations of observatories. 

In any case, inserting a ~ 2 ir£ H /£ B ~ 10 3 ” 4 and T max ~ 10 6 GeV into eqs. (6.9), (6.10), 
the thermal background would be ~ 40 orders of magnitude below the desired level at the 
peak eLISA frequency (cf. also fig. 2). What is significant about the thermal background, 
though, is that it continues to grow with k for another more than 10 decades, and therefore 
eventually overtakes the decaying non-equilibrium signal at short distance scales. In fact, 
originating as it does from fluctuations at the scale k ~ 3T max and red-shifting as dictated 
by entropy conservation, the peak power is in the range k 3T dec (3.9/106.75) 1 / 3 T dec at 
the time of photon decoupling, and in the microwave range today. Therefore it falls in the 
range of recently conceived high-frequency experiments [53]—[57]. 


8. Conclusions and outlook 

We have estimated the magnitude and shape of the gravitational wave background that 
is produced by Standard Model physics during the thermal history of the universe until the 
temperature To ~ 160 GeV corresponding to the electroweak crossover, cf. eqs. (6.2), (6.8) and 
fig. 2. The infrared part could have been of potential interest in that the forthcoming eLISA 
experiment is probing sub-Hz frequencies with unprecedented precision. Unfortunately, we 
have found that in this range the thermally produced gravitational wave signal is many orders 
of magnitude below the observable level, cf. sec. 7. 

In general, the thermally produced gravitational wave background resembles a bit the 
blackbody spectrum of photons and neutrinos. Its shape is not the same because gravita¬ 
tional waves never equilibrate. Therefore the shape has to be determined by a dynamical 
computation which has not been carried out even at full leading order. Nevertheless, it is 
already clear from a leading-logarithmic estimate that the peak of the power today lies in 
the same microwave domain as for photons and neutrinos (cf. fig. 2). Therefore, the best 
observational prospect lies with high-frequency (0.1 - 4.5 GHz) experiments [53]—[57]. At the 
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current stage it seems challenging to reach a sensitivity below fl GW ~ 10 5 [56], whereas an 
optimistic theoretical expectation would be 


0 0 1 ( 4.5 GHz 

^GW~°radX 10Q X ^ 10() GHz 


3 


10 


-li 


( 8 . 1 ) 


where 1/100 corresponds to the maximally allowed fraction in gravitational waves at the 
electroweak epoch when all degrees of freedom are relativistic, and 100 GHz to the frequency 
associated with generic blackbody radiation. So, there is surely a long way to go till detection. 

There is, however, one consideration which can already be carried out. Indeed, unlike 
neutrinos, gravitational waves must not carry as much energy density as one relativistic 
degree of freedom [51, 52]. This constrains the total energy density stored in them and, 
given that the production rate peaks at the maximal temperature, the maximal temperature 
reached. The total energy density corresponding to eq. (6.8) can be estimated as 


J dln&o fI GW (&o) - 


24?) T max 

7 tV 67re 3 m P1^0 
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_ 
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V \ ^max 
3 ) m pl 
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where we varied <f> between two limits: the factor 8 originates if we adopt the form of </> 
appearing on the second line of eq. (6.2), setting the unknown constant to zero and the 
couplings to values mentioned in the caption of fig. 2, whereas the factor ?)/3 originates if we 
use the first line of eq. (6.2) and cut off the integral at ko = Tq: J q T ° dko fcg = Tq /3. According 
to Planck data [58] only a small fraction of a relativistic degree of freedom beyond those in 
the Standard Model can be permitted, so at Tq ~ 160 GeV we must require 


24 /g /)\ Tmax ^ 1 

7 tV 6ire 3 V 3 J 77? pi 100 

Inserting e ~ 35, ?) ~ 400 we obtain T max < 10 17 ...10 1S GeV. This is not a very strong 
constraint, 3 but the estimate could be sharpened with more knowledge about the function cj). 

To summarize, a determination of the function cj), defined by eq. (6.1), beyond the leading- 
logarithmic terms that we have obtained here, seems to pose an interesting problem. This 
computation represents a well-defined challenge in thermal field theory, analogous to that 
for the photon production rate from a QCD plasma [39, 40] or the right-handed neutrino 
production rate from a Standard Model plasma [43, 59]. It is technically more challenging, 
because every single particle species carries energy and momentum, and therefore we leave 
the practical implementation to future works. (A somewhat related computation, but for the 
off-shell kinematics k = 0, uj>T, has been presented in ref. [60].) 

3 In particular, within the standard inflationary paradigm, reheating temperatures above ~ 10 16 GeV are 
considered all but excluded. 
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